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Abstract 

Three phenomena are involved in sand movement: erosion, wind transport, and sedimentation. This paper presents a 
comprehensive easy-to-use multiphase model that include all three aspects with a particular attention to situations in 
which erosion due to wind shear and sedimentation due to gravity are not in equilibrium. The interest is related to the 
fact that these are the situations leading to a change of profile of the sand bed. 
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1. Introduction 

When the shear stress exerted by wind on a sandy surface is sufhcietly strong, sand grains are lifted from the 
sand bed and are transported by wind to sediment downstream. The raising sand grains follow a ballistic trajectory 
influenced by drag and gravity, eventually impacting again on the surface and inducing new particles to detach from 
the surface. This phenomenon, knows as saltation, generates a layer close to the sand bed with a typical maximum 
height of 10-20 cm. Saltation is the main reason of erosion of sandy surfaces and together with the consequent 
sedimentation of sand particles it is the main reason of dune motion and accumulation of sand in specific regions 
where recirculation occurs. 

The engineering interest in understanding and simulating the dynamics of windblown sand, e.g. dune fields of 
loose sand, is dictated by their interaction with a number of human infrastructures in arid environments, such as roads 
and railways, pipelines, industrial facilities, farmlands, towns and buildings as shown in Fig. up. 

Moving intruder sand dunes, soil erosion and/or sand contamination can be comprehensively ascribed, from a 
phenomenological point of view, to non-equilibrium conditions, where the two processes, erosion and sedimentation, 
do not balance, leading to the erosion or deposition of sand on the soil and eventually to the evolution of that interface. 
In other terms, such non-equilibrium situations are the most interesting cases from the applicative point of view. 

Bagnold (8) was the first who studied sand erosion and postulated a relation for the sand flux, determining the 
importance of wind speed and of the related shear stress on the sand surface. Later authors ll5l fl2l [T7112T1 [26l [ 281 [30l (31 3 
introduced several corrections to Bagnold’s rule, but all the models have in common the observation that a sand grain 
is ejected from a sand bed if and only if the shear stress at the surface is larger than a threshold value. 

Sauermann et al. l30l observed that saltation reaches a steady state after a transitory phase of 2 seconds. After 
this period the trajectories are statistically equivalent for the ensemble of grains. This phenomenon happens because 
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Figure 1: Windblown sand interaction with anthropic activities: farmlands (a), towns (b), railways (c) 


the new ejected particles increase the sand concentration in the saltation layer and this reduces the speed of saltating 
grains. So, a steady state is reached when all particles are ejected with the same velocity (see also E1I30I3D). A 
nice mathematical models of the saltation phenomenon is proposed by Herrmann and Sauermann m who studied 
the dynamics of the surface of a dry granular bed dividing the sand bed into a non-moving time-dependent region 
providing sand mass and another time-dependent region above it in which sand particles can move transported by the 
wind. They propose a model averaged over the vertical coordinate, presenting a free boundary. 

Ji et al. H3 coupled a k — e model with a multiphase approach in which the slip of the dispersed phase is modeled 
by an algebraic model. Similar turbulent one-dimensional models are proposed in B24l [27 ;. 29 301, however without 
a multiphase coupling. Kang and coworkers IfTSl fl9] i20l instead couple a multiphase model for the fluid flow with a 
particle method for the sand grains. A similar coupling was also used in Il2ll4l l23l where however the wind flow was 
computed indipendently from the presence of sand particles via a suitable turbulence model, typically the k-e model. 

Sedimentation has also been widely studied in the literature starting from several applications mainly in envi¬ 
ronmental and chemical engineering. One of the most important component in this phenomenon is the drag force 
experienced by the sedimenting particles that has driven a lot of attention by many authors as well reviewed in |9fl . 

Differently from previous papers, here we will propose a comprehensive multiphase model for the entire process 
including sand erosion, wind transport, and sedimentation, that working also in non-equilibrium conditions is able to 
deal with the development of the stationary saltation layer starting from generic initial and boundary conditions and in 
particular from clear air and oversaturated situations. In order to do that we develop a so-called first order model (in 
time) of sand erosion, transport and deposition, that can be easily tuned using experimental test cases. The resulting 
advection-diffusion equation for the suspended phase can then be coupled with a k- a> model describing the turbulent 
fluid flow. The mathematical model can then be solved with the aid of the fundamental erosion/deposition boundary 
condition at the sand bed, that depends on the shear stress. 

The plan of the paper is then the following. After this introduction. Section 2 presents the mathematical model 
mainly focusing on the advective phenomena, on the microscopic dynamics related to the collision between sand 
grains, and on the erosion boundary condition. The result of some numerical simulations focusing on how the station¬ 
ary condition is reached when wind blows over a heterogeneous sand bed are reported in Section 3. 

2. The Erosion/Transport/Deposition Model 

We consider the flow of sand as a multiphase system composed of sand grains in air. Single sand grains have a 
density p s and float in air with a volume ratio <f> s (typically well below 1%), so that the partial density of sand in air is 
p s = p s (p s . Saturation obviously implies that <pf = 1 - <f> s where 0/ is the volume ratio of air. The mixture of air and 
sand grains is flowing on a sandy surface having a close packing volume ratio <j> s . 

Because wind flow is in a turbulent regime the fluid phase is modelled by the Reynolds-averaged Navier-Stokes 
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equations (RANS) equations. More precisely, a k - to turbulence model is selected to provide the closure l25l 
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with standard boundary conditions. In ([T]) A: is the turbulent kinetic energy, to is the specific dissipation rate, v„ and v, 
are, respectively, air and turbulence viscosities, P k and P,„ are the production terms for k and to, and y and C w are two 
empirical costants. 

In describing the transport of sand we start observing that while sand particles are trasported by the wind they drift 
down with a characteristic sedimentation velocity due to the action of gravity 0. In addition, particle collide giving 
rise to an extra-flux term q co //. Hence, one we can write the following equation for the sand volume ratio 
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This closure can be actually deduced under suitable modelling assumptions from a more general multiphase model 
involving mass and momentum balance for the suspended phase. 

The sedimentation velocity w se d can be evaluated by the balance of drag and buoyancy forces and strongly depends 
on the grain size. For instance, if we define the particle Reynolds number as the one felt by the sand grains of diameter 
d during their flow and therefore based on the relative velocity between air and solid particles, Re s = <t>f\vf - v s \d/vf, 
then in the so-called Newton regime, corresponding to particle Reynolds numbers above 500, the drag coefficient Cd 
is approximately constant (for instance, Cd = 3 is used in |[27l[30i ). so that one has the classical relation 
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However, at the other extreme, i.e., for particle Reynolds number below few units, corresponding to the so-called 
Stokes regime, Cd ~ jy ■ so that one has the classical Stokes sedimentation velocity 


W sed 
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Considering that in aeolian sand trasport the phenomenon is limited to the first few centimeters from the ground and 
that there the particle Reynolds number is below 100 (see Fig. [2^), a better evaluation of the sedimentation velocity 
with respect to Eqs. ([4] [5]) can be obtained fitting the experimental dependence of the drag coefficient on the particle 
Reynolds number as reviewed in 0 and shown in Fig. [ 2 J 5 . 

Coming to the collision term q eo //, already introduced by Batchelor Da, neclecting it would imply that sand grains 
are only transported under the action of drag and gravity. However, collisions among particles have the important 
non-negligible effect of generating a sort of diffusion of sand particles from higher to lower density areas, that results 
fundamental in this modelling framework for the stationary formation of the saltation layer. 

For high volume ratios near close packing, Auzerais et al. f7[ suggested the following nonlinear law on the basis 
of experimental data 

q coil = -D eff V= -— , with k e [2,5]. (6) 
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Figure 2: (a) Particle Reynolds number as a function of the sand grain diameter of a sedimenting particle in air. (b) Sedimentation velocity of sand 
grains in still air in the Stokes limit of Eq. (5) and according to the experiments summarize in {9). (c) Qualitative dependence of the diffusivity 
coefficient v e yy from the distance from the sand bed. 


Such a term enforces the need of avoiding that the close packing volume ratio cp s is reached for the sedimenting mass. 
However, as <p s <s <p s in wind-blown sand applications, the relation can be simplified to 

Qco// = ~^effps • (7) 

From the constitutive viewpoint, the collision terms 0 or 0 can be considered as deriving from treating the ensemble 
of particles as a gas, so that the stress term for the solid constituent is isotropic through a coefficient that depends on 
the particle density. Substituting Eqs. ([7]) and ([3]) into (|2ji gives the nonlinear degenerate advection-diffusion equation 

^ + V ■ 0 p s \ f ) - ^-(<p,w sed ) = V • {y eff <p k s ^W(ps) . (8) 

Actually, if the limit value k = 1 is also allowed, one has the linear case sometimes used in the literature. 

The coefficient v,,ff might be considered as composed of three contributions v e ff - v s + v mr t, + v„,„;, that take into 
account of the possible dependence from 

• the shear rate, or better the velocity gradient v s = Vj(V0,); 

• the turbulence of the flow, so that this term results from the integration of the CFD simulation in a turbulent 
regime; 

• the molecular diffusion, but as reported in ffT8l , this term can be neglected with respect to other quantities. 

Actually, starting from the obvious observation that the behaviour of sand particles is isotropic. Objectivity, i.e., 
independence of the constitutive dependence from the reference frame, implies that v s is a scalar isotropic function 
of a tensor. By the representation theorem of isotropic function v s can then only depend on the invariants of the rate 
of strain tensor D = yiVv + Vv r ). However, since the flow can be considered as a perturbation of a shear flow in 
the vertical plane, the leading contribution is the second invariant Ed — y |(trB) 2 - trO 2 ]. For this reason we assume 
that v s = v s (Ed). This is an important generalization because all papers refers to a dependence on the wind velocity 
u*, which is a well defined quantity only on horizontal surfaces, while more general surfaces like dunes require a 
dependence from an objective invariant of the rate of strain tensor. 

In order to understand the meaning of the collision term we can look for stationary configurations for which all 
quantities depend only on the quote z and velocities are directed along a flat plane (at z = 0 in the direction x). In this 
case, the stationary profile in the saltation layer can be obtained integrating 

, _, dip s 

v eff<P s + W ^ed(ps = 0 . 

jointly with the boundary condition at the sand bed. In the simplest case in which v c yy is constant and k = 1, one 
immediately has an exponential profile with a characteristic length related to the thickness of the saltation layer given 
by 6 S = v e ff/ws e d and an integration constant related to the erosion boundary condition. 

In general, referring to Fig. [2J; the dependence of the coefficient on the distance from the sand bed shows a strong 
diffusion closer to the surface and a low diffusion at a distance close to the height of the saltation layer, while particles 
that escape from the saltation layer are again strongly mixed due to the increasing diffusion due to turbulence. 
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The general features of the erosion boundary condition can be inferred from experiments known for a long time 
that show that, generally speaking, erosion only occurs if the shear stress at the surface r exceeds a threshold level r t , 
or equivalently that u* = sjr/pf exceeds a threshold level u* = jT,/pf (see li22l and referencces therein). 

Referring to the last notation, because it is the one classically used in the literature (though from the numerical 
point of view what is computed is t which is then compared to r,) one has the flux boundary condition 

<7s = -v e/ /0j“'V0 v • n = p(u*~ - , 

where (/) + = stands for the positive part of /. Bagnold’s formula [8J u* = Cf ’ gd is quite successful in 

determining the dependence of u* from the grain diameter. On the other hand, the quantification of the parameter /3 
is not straightforward because most of the experiments measures the horizontal sand flux parallel to the surface while 
for the boundary condition one would need some knowledge of the vertical flux perpendicular to the surface which is 
the one related to the ejection of sand grains from the sand bed. Very recently, on the basis of their experiments Ho et 
al. 114] 15 16] proposed /3 = AhP/ sjd/g, so that the sand flux due to erosion takes the form 


= 


wapf 

gd 




(9) 


where w = 0.5 m/s is the sand grain ejection vertical velocity evaluated experimentally [1.1, 14j and a is a dimension¬ 
less free parameter to be fitted to experimental sand flux profiles. 


3. Simulations 

As domain of integration we focus on the flow over a horizontal heterogeneous lane. This is neither an artifi¬ 
cial situation, nor a case of limited importance. In fact, most of the landforms in arid regions and roads are well 
approximated by a horizontal flat plane. This geometrical setup is also retained in most of the wind tunnel experi¬ 
mental studies present in the literature that are however mainly addressed to the characterisation of the windblown 
sand concentration and flux in uniform, in-equilibrium steady state conditions. Nevertheless, uniform and steady state 
conditions are excessively ideal ones. In these situations the incoming wind already transports the maximum allow¬ 
able sand density, so that erosion and deposition are balanced, and hence the sand bed surface is neither scoured, nor 
accumulated. Conversely, in many engineering applications, attention must be paid to non-equilibrium conditions that 
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Figure 3: Considered setup with alternation of erodible and non-erodible surfaces, of equilibrium and non-equilibrium regions, and a sketch of the 
height of the saltation layer 6 S . 


will cause erosion or settlement of the sand bed. So, as sketched in Fig. [3] our horizontal plane is characterized by the 
alternation of erodible sandy regions and non-erodible regions, corresponding for instance to a street. In the simula¬ 
tion the inflow wind is clean and with a fully developed logarithmically shaped velocity profile, with zo = ICC 5 m and 
u* ranging between 0.3 and 1 m/s. The threshold value for erosion is u* = 0.25 m/s and the grain size is d = 0.25 mm. 

From the simulation shown in Figures[4]and[5] as soon as wind overcomes the boundary between the non-erodible 
and erodible surface (that is put at a distance L w = 1 m from the inflow boundary) the saltation layer starts to develop. 
Figure [4ji plots the profiles of the sand volume ratio <p s at several points at the beginning of the erodible zone. The 
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Figure 4: Sand density profiles in the saturation (a) and in the sedimentation (b) regions 
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model correctly predicts the progressive uptake of sand and increase in the depth of the saltation layer, till saturation 
is reached because of the equilibrium between erosion and sedimentation. The thickness of the saltation layer at 
equilibrium is about 10 cm in qualitative agreement with experiments. Viceversa, as shown in Fig. |4j), at the beginning 
of the second non-erodible zone, there is a reduction of the windblown sand density in the upper part of the stream 
because of the sedimentation process, not balanced by saltation. 

Actually, due to diffusion, some sand also diffuses upstream, mainly very close to the surface where diffusion 
is dominant. Referring to Fig. |5};, in this region one can notice an exponential growth of the scaled total sand flux 
Q! Qsat where Q sa , is the sand flux at equilibrium. After the first soil discontinuity at x = 0 convection dominates 
and Q saturates in a length close to L sat = 3 m (see Fig. [5j)), that is nearly independent from the scaled wind velocity 
u* /u*, in qualitative agreement with |6j . 

When the erodible surface ends the total sand flux Q readily decreases in a characteristic distance that increases 
with the wind velocity as shown in Fig. [5j>. It can be noticed from Fig. [5jl that the decrease is less than exponential in 
qualitative agreement with tm 

In conclusion, the model m jointly with the fundamental erosion boundary condition ([9]> and other standard 
initial and boundary conditions is able not only to describe the erosion/transport/sedimentation process in stationary 
situation, but also to capture all features characterizing the development of the saltation layer up to equilibrium and 
of its reduction, that occur in nature when the sandy surfaces is heterogeneous. 
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Figure 5: Scaled total sand flux over the sand bed at different wind velocities 
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